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Introduction 

Steady flow over the leading portion of a multicomponent airfoil section is studied us- 
ing computational fluid dynamics (CFD) employing an unstructured grid. To simplify the 
problem, only the inviscid terms are retained from the Reynolds- averaged Navier-Stokes 
equations — leaving the Euler equations. The algorithm is derived using the finite-volume 
approach, incorporating explicit time-marching of the unsteady Euler equations to a time- 
asymptotic, steady-state solution. The inviscid fluxes are obtained through either of two 
approximate Riemann solvers: Roe’s flux difference splitting or van Leer’s flux vector split- 
ting. Results are presented which contrast the solutions given by the two flux functions as 
a function of Mach number and grid resolution. Additional information is presented con- 
* Research Engineer, Aerothermodynamics Branch, Space Systems Division. 
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cerning code verification techniques, flow recirculation regions, convergence histories, and 
computational resources. 


Nomenclature 

A ceil Cell area 
ei/t Wave speed 
c Speed of sound 
e Specific energy 
F, G Cartesian inviscid flux vectors 
h Specific enthalpy 
hk Component of H vector 
H Locally face-normal flux vector 
Lqo Infinity error norm 
L2 Least-squares error norm 
M Mach number 
p Pressure 

R Matrix of wave signature column vectors 
Rk Column vector component of R 
As Length of cell side 
t Time 
At Time step 

T Transformation matrix from Cartesian coordinates 
u,v Cartesian velocity components 
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U Conservative variables 


AV Wave strength vector of jumps 
x,y Cartesian coordinates 

to local normal/tangential coordinates 

a m Runge-Kutta coefficient 
7 Ratio of specific heats 

6 Polar angle 

v Courant number 

p Density 

<I> Flux function 


Subscripts 

( ) 0 Image cell value 
( Interior cell value 

( ) k Vector component (k=l,2,3,4) 
( ) L States in the “left” cell 
( States in the “right” cell 

( ) Normal component 
( ) t Tangent component 
( j^Freestream quantity 
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Superscripts 

( )” > ( )" +1 Time levels 

( * ) Roe-averaged quantity 

( )* Smoothed values 

( ) + Right-traveling information 

( )~ Left-traveling information 


Governing Equations 


Two-dimensional, unsteady flow is described by the following conservation laws 




where the state vector is given by 


U 


/ \ 

P 

pu 
pv 

\ pe ) 
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and the inviscid fluxes are 


f \ f \ 


pu 


pv 

pu 2 -fi P 

G = 

puv 

puv 


pv 2 + p 
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The equation of state closes the system 


e = 


P 


(7-1 )P 2 


+ \ (« 2 + v 2 ) 


and for convenience, enthalpy is defined as 


h = 



Note: for air: 7 = 1.4. 


Flux Functions 

The inviscid fluxes were computed using two different upwind schemes: flux difference 
splitting (FDS) of Roe 1 and flux vector splitting (FVS) of van Leer . 2 Both techniques are 
outlined below. 
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van Leer’s Flux Vector Splitting 

For van Leer’s flux vector splitting scheme, the inviscid portion of the fluxes across a 
cell face is given by 

* (U L , Ur) = T- 1 (H+ (U l ) + H- (Ur)) 

where T -1 is the inverse transformation from cell face oriented coordinates to Cartesian 
coordinates 

10 0 0 
0 sinO cosQ 0 

T 1 = 

0 —cosO sinO 0 
0 0 0 1 

The fluxes are split into right-traveling and left-traveling based on the cell-centered, 
face-normal Mach number, i.e., 

M n > 1 H + = H (U L ) H~ = 0 

M n < -1 H + = 0 H- = H(Ur) 

\M n \ < 1 H+ = H+(U l ) H- = H-(Ur) 

where 

pu n 

pul+p 

pu n u t 
pu n h 
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and 



where the wave speeds are 



the jumps in the wave strengths are 
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and the wave signatures are given by 


R — R-3 Ra ^ 

/ . . \ 


1 

0 

i 

1 

ii — c sinO 

c cosO 

u 

u + c sinO 

v + c cosO 

c sin6 

V 

v — c cosO 

h — u n c 

UtC 

l (« 2 + V 2 ) 

h + u n c 


where A ( ) represents the jump between the left and right states 

M ) = ( )«-< )u 

and the ( ' ) quantities are the Roe-averaged variables 


P = y/pLpR 


it = 


v — 


h = 


y/pZ^L + y/pRUR 

S/Pi + y/pR 
y/PLVL + y/pRVR 
y/pL + y/pR 
y/pZ^L + y/pRh^R 

yffZ + yfpZt 


where c, u n , and u t are calculated directly from p , u, u, and h, so 


p = iizStl 
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u n = usinO — vcosO 
u t = iicosO + vsinO 


To prevent expansion shocks, an entropy fix is imposed. A smoothed value of |a^| is 
defined for the acoustic waves ( k = 1 and k — 4 ) 


H* = 


|a*| \a k \ > \6a k 

ik + \ Sa * M < ¥ a * 


with 


8c ik = max(4Aa^,0) 


This provides a parabolic (and thus continuous) curve where the wave speeds change signs 
(e.#., in a transonic expansion, or at a stagnation point). 


Time Integration 

Time integration of the governing equations was performed by two methods: forward 
Euler and Runge-Kutta, with the time step per cell area was computed via 


At v 

A ce ll max faces [(u n + c) As] 


where v is the Courant number. 
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Forward Euler 


The simplest scheme is that of forward Euler 


U n+I = U” 


At 

A cel I 


R(U n ) 


where R is the residual of the cell given by 


R(U n ) = * As 

faces 

This scheme was only used as a step in the debugging process enroute to the following 
multi-stage scheme. 

Runge-Kutta 

A four-stage, optimally-smoothing, Runge-Kutta, time-stepping scheme due to Tai 3 was 
implemented as follows 


U° 

U 1 

u 2 

u 3 

u 4 


u n+1 


u n 

U °-“'£; R K> 

u ° 

U°-o 3 ^-R(u 2 ) 

Acell V ' 

u»- q< ^-r(u 3 ) 

Acell V ' 

U 4 


10 


with a k = [0.0833,0.2069,0.4265,1.0]. The scheme was nominally run with a Courant num- 


ber of 2.0. 


Boundary Conditions 

The solid wall boundary condition for the airfoil was enforced in a weak sense by setting 
the cell-centered state in an image-cell located just inside the solid boundary surface. Physi- 
cally, flow tangency must be preserved at the wall, as well as a zero pressure gradient normal 
to the wall. This image-cell wall-boundary procedure is of first order accuracy, because it 
neglects the wall curvature. This will introduce inaccuracy in the case of a highly curved 
wall that is not resolved by sufficiently fine cells on the wall. This defines the values in the 
image-cell as follows 


P 0 

= Pi 

Po 

= Pi 

( u t)o 

= K), 

( U n)o 

— i u n) 


where the ( ) 0 represents an image-cell quantity and ( is the appropriate interior cell 
quantity. The flow tangency condition is supplied by merely reflecting the normal component 
of velocity across the boundary face. In more realistic terms, this boundary condition can 
be thought of as a symmetry-plane condition. 

For the farfield boundary, the image-cells were specified as freestream conditions, and 
the Riemann solver simply picks the proper information to use depending on whether the 
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local flow is into or out-of the computational domain. Note: a more careful outer boundary 
condition for a lifting airfoil would superimpose a potential vortex velocity field at the outer 
boundary. 


Computational Meshes 

Three different sizes of meshes were used for this project as indicated in Table 1. The 
results were computed on the medium grid unless otherwise specified. Figure 1 shows a view 
of each of the grids in the vicinity of the airfoil. 

The Debugging Process 

Debugging the code is the most labor intensive part of any project — if it is not done in 
a logical manner. .Steps along the debugging path were as follows: 

1. Check pointers for a simple mesh. 

2. Verifying cell areas were positive and correct magnitude. 

3. Implement only forward Euler time stepping. 

4. Setting all boundary conditions to be freestream conditions. 

5. Verify both flux functions agree with one another. 

6. Run angle of attacks: 0°, ±45°, ±90°. 

7. Add solid wall boundary conditions. 

8. Run angle of attacks: 0°, ±45°, ±90°. 
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9. Turn on the multistage time stepping. 


10. Run angle of attacks: 0°, ±45°, ±90°. 

11. A freestream Mach of 2.5 was run to determine (very readily) that the flow was going 
in the right direction-indicated by a very pronounced bow shock encompassing the 
front of the body. 

12. A grid convergence study was conducted to check for consistency of the algorithm-see 
following section. 

The final check was done just from knowledge of how the flow should behave-common sense. 

Results and Discussion 

The flow over the leading section of a multicomponent airfoil is computed for three 
different freestream Mach numbers: 0.4, 0.8, and 1.2. The Mach 0.8 case was arbitrarily 
chosen for a grid convergence study, and it was run on all three meshes: coarse, medium, 
and fine. Each Mach number/grid combination was also run using both of the flux functions 
discussed previously. 

Computational Resources 

The code was primarily run on a Hewlett Packard Apollo series 700 workstation. When 
compiler-optimized, the code typically ran around 2.6 x 10“ 4 CPUs/iteration/cell when using 
the FDS flux function. For example, a 6124 cell mesh running for 2200 iterations would have 
a total CPU time of around 1 hour. Timings for different machines and flux functions are 
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given in Table 2. As coded, the FVS scheme runs 40 percent faster than the FDS scheme. 
Note that the Sun Sparc station runs an order of magnitude slower than either of the HPs. 

Mach Number Effects 

Figure 2 depicts Mach contours of the flow in the vicinity of the airfoil for three Mach 
numbers: 0.4, 0.8, and 1.2. The sonic line and the M = 0.5 contour lines are labeled in each 
subfigure. The remaining contour lines occur at intervals of 0.1 Mach. The results in the 
left column of the figure were produced using Roe’s FDS flux function and those in the right 
column resulted from van Leer’s FVS flux function. Comparing the left and right columns 
of Fig. 2 it is apparent that the two different flux functions give nearly identical results 
for the flow field. The only differences are minor and occur in the trailing edge/flap cavity 
region which will be investigated further in the following section. The Mach 0.4 case shows a 
stagnation region at the leading edge followed by two expansions: one near the leading edge 
on the suction side of the airfoil and another just before the flap slot on the pressure side. 
The Mach 0.8 flow contains two shocks: a strong one on the upper surface around 90 percent 
chord and a smaller one following a transonic expansion in the flap cavity. There is also a 
significant supersonic “bubble” on the upper surface. At Mach 1.2, the upper surface shock 
moves to the trailing edge and forms the familiar “fish-tail” shock structure. Also apparent 
is a weak, detached bow shock standing well off the leading edge. 

Flap Cavity Region 

The differences between the two flux functions become more apparent upon closer in- 
spection of the flap cavity region. Figure 3 shows streamlines around the trailing portion of 
the airfoil for a freestream Mach of 0.8 on all three grids for both flux functions. The figure 
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shows the results using Roe’s FDS in the left column and results using van Leer’s FVS in 
the right column. Notice that for all three grids, FDS shows a vortical structure in the flap 
region, while FVS only shows a vortical structure for the two finer grids. In addition, FDS 
supports a much more complex structure — the fine mesh showing three interacting recircu- 
lation regions. This is a result of the lower dissipation inherent in Roe’s FDS scheme as 
compared to van Leer’s FVS scheme. 

Grid Convergence 

Shown in Fig. 4 are Mach number contours about the airfoil for three different grids. The 
flow has a freestream Mach number of 0.8 and was computed using the FDS flux function. 
As portrayed in the Fig. 4, the global solution only changes with respect to the thickness of 
the shock, implying that the scheme is consistent with respect to the governing equations. 

Convergence Histories 

A summary of the total number of iterations required to reach an L 2 error norm of 
the energy equation of 1.0 x 10 -6 is shown in Table 3. Roe’s FDS scheme takes longer to 
converge than van Leer’s FVS scheme in every case. This is apparently due to the highly 
dissipative nature of FVS which tends to smooth spurious transients. In two cases, the FDS 
scheme even fails to converge to the specified error tolerance-more on this to follow. 

Shown in Fig. 5 are normalized convergence histories for three different Mach numbers 
using the two flux functions. (Note: the straight line at the tail of the iteration line-plots is 
an artifact of the plotting routine used.) This figure clearly shows the convergence problem 
inherent in this application of Roe’s FDS flux function. Notice the cyclic behavior of the 
error residual. This corresponds to the following cycle: a flow feature moving just slightly, 
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the disturbance propagating into the rest of the field, and a reflected disturbance telling 
the feature to move back. The process repeats with a cycle corresponding to speed of the 
disturbance and the number of cells through which it passes. This was verified by halving 
the Courant number and noting a increase by a factor of two in the cycle period — the shape 
remained the same. 

Shown in Fig. 6 and Fig. 7 are the regions which contained the most error norms 
during the run for each of two cases. The left side of each figure is a blocked-contour plot of 
the number of times a particular cell was responsible for the error norm. The right side of 
the figures shows the streamlines in the same viewing area. Figure 6 corresponds to a Mach 
0.8 flow using the fine grid and Fig. 7 corresponds to a Mach 1.2 flow using the medium grid. 
For the Mach 0.8 case, the flow features responsible for the convergence problem appears to 
be both the transonic expansion on the bottom edge of the flap cavity region and the shock 
standing on the upper surface of the airfoil. However, in the Mach 1.2 case, the flow feature 
responsible for the convergence problem is the shock just ahead of the recirculation zone in 
the center of the airfoil flap cavity region. 

An attempt was made to alleviate the convergence problem by “smoothing” the grid 
slightly. This was done by telling each node to move toward the centroid formed by its 
neighbors. This smoothing process altered the grid enough to stabilize the flow features in 
slightly different locations — hopefully allowing convergence; but the residual just hung at a 
slightly lower error norm. 
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Conclusion 


The results indicate that an airfoil experiences radically different flow fields as a function 
of Mach number. Features range from smooth variations of the flow properties at a low Mach 
numbers to discontinuous shocks forming at higher Mach numbers. In spite of the range of 
flow conditions, both flux functions appear to work quite well. Even though FDS is slightly 
more expensive than FVS and it experiences convergence difficulties, it appears to do a better 
job resolving recirculation regions than does FVS. 
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Table 1: Mesh statistics 


Name 

Cells 

Edges 

Nodes 

coarse 

2697 

4098 

1401 

medium 

4025 

6114 

2089 

fine 

6124 

9308 

3184 


18 




Table 2: CPU seconds per iteration per cell 


Machine FDS 

Sun Sparc 1+ 3.9 x 10 -3 

HP Apollo 710 2.6 x 10~ 4 
HP Apollo 720 2.5 x 10~ 4 


FVS 

1.4 x 10~ 3 
1.8 x 10~ 4 
1.7 x 10~ 4 



Table 3: Number of iterations for convergence 


Mach 

Grid 

FDS 

FVS 

0.4 

medium 

2542 

2123 

0.8 

coarse 

1738 

1496 


medium 

2733 

1909 


fine 

hung 

2396 

1.2 

medium 

hung 

1823 




(c) Fine grid. 


Figure 1: Computational meshes employed. 
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(c) Mach 0.8 FDS. (d) Mach 0.8 FVS. 



(e) Mach 1.2 FDS. (f) Mach 1.2 FVS. 


Figure 2: Comparison of Mach number contours for three Mach numbers (AM = 0.1) 























Figure 6: Regions responsible for the convergence difficulties of FDS scheme at Mach 0.8, 

fine grid. 
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(a) Blocked-contours of Loo norm occur- (b) Streamlines in the same vicinity, 

rences. 


Figure 7: Regions responsible for the convergence difficulties of FDS scheme at Mach 1.2, 

medium grid. 
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